% This code simulates the model for EIS ~= 1. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SIMDATA = Sim_EIS(l_q , l_q_burn , N_sim , nu , phi , psi , c_0 , Var_Names)

% INPUT:
% -- l_q: Number of sample periods
% -- l_q_burn: Number of burn-in periods
% -- N_sim: Number of simulations
% -- nu: constant gain parameter
% -- phi: informativeness of prior
% -- psi: EIS parameter
% -- c_0: starting value of log real per capita consumption level
% -- Var_Names: names of variables to be stored
%
% OUTPUT:
% -- SIMDATA: stored simulations as a cell

%% 1. LOAD PARAMETERS
if nu == 0
    error('nu cannot be zero!')
end

[I , J] = param_EIS(nu , phi , psi);

%% 2. SIMULATIONS
% additional parameters for simulation
T = l_q;                %Number of sample periods
T_burn = l_q_burn;            %Number of burn-in periods
T_total = T_burn + T;     %Total periods including burn-in periods for simulation
d_0 = c_0 + I.d_c_mean;   %Starting value of log dividend

%% 2.1 Simulate the economy
% Simulate quarterly iid consumption growth, column for each simulation
cg_full = normrnd(I.mu , I.sigma , T_total , N_sim); 

% Simulate quarterly shocks to dividend growth
dg_shock = normrnd(0 , I.sigma_d , T_total , N_sim); 

% Log real consumption level from 1 to T_total
c_full = c_0 + cumsum(cg_full); 

% Log real dividend level from 1 to T_total
d_full = zeros(T_total , N_sim); 

d_full(1 , :) = I.lambda * cg_full(1 , :) + (1 - I.alpha) * d_0 + I.alpha * c_0 + I.alpha * I.mu_dc + dg_shock(1 , :);

for k = 2 : T_total
  d_full(k , :) = I.lambda * cg_full(k , :) + (1 - I.alpha) * d_full(k - 1 , :) + I.alpha * c_full(k - 1 , :) + I.alpha * I.mu_dc + dg_shock(k , :);
end

% Log dividend growth from 1 to T_total
dg_full = [d_full(1 , :) - d_0 ; d_full(2 : end , :) - d_full(1 : end - 1 , :)]; 

% Four-quarter change
cg_yoy_full = (c_full(5 : end , :) - c_full(1 : end - 4 , :))/4;
dg_yoy_full = (d_full(5 : end , :) - d_full(1 : end - 4 , :))/4;

% Construct State Variable \tilde{\mu}
tmu_full = expweightedavg(cg_full , I.nu);

% Calculate risk-free rate
rf_full = -I.tmu_m_s + tmu_full * I.phi/I.psi - 0.5 * I.xi_s^2 * I.sigma^2;

%Rf_full = exp(rf_full) - 1; 

%% 2.2 Calculate price-dividend ratio and price
pd_full = J.Ad_0 + J.Ad_1 * I.phi * tmu_full + J.Ad_2 * (d_full - c_full);

P_full = exp(pd_full + d_full);

% Log return implied by log-linearization
ret_ll_full = J.kappad_0 ...
            + J.kappad_1 * pd_full(2 : end , :) - pd_full(1 : end - 1 , :)...
            + dg_full(2 : end , :);
        
% Log return calculated by price and dividend
ret_full = log( P_full(2 : end , :) + exp(d_full(2 : end , :)) ) ... % log(P_{t+1} + D_{t+1})
         - log( P_full(1 : end - 1 , :) ); % - log(P_t)
     
ret_MS_full = movsum(ret_full , 4 , 'Endpoints' , 'discard');     
     
% Excess returns     
EXRET_full = exp(ret_full) - exp(rf_full(1 : end - 1 , :)); % Net excess returns
EXRET_ll_full = exp(ret_ll_full) - exp(rf_full(1 : end - 1 , :)); % Net excess returns

exret_full = ret_full - rf_full(1:end-1,:); %log excess returns
exret_ll_full = ret_ll_full - rf_full(1:end-1,:); %log excess returns

% Subjective expected log excess returns    
sub_ret_full = I.phi/I.psi * tmu_full + (J.kappad_1 - 1) * J.Ad_0 + J.kappad_0 + (J.kappad_1 * J.Ad_2 + 1) * I.alpha * I.mu_dc ...
             + (1 - I.phi) * I.mu * (J.kappad_1 * J.Ad_1 * I.phi * I.nu + (I.lambda - 1) * J.kappad_1 * J.Ad_2 + I.lambda); 

% Experience returns
tmur_1_full = expweightedavg(ret_full , I.nu);   %2 to T_total, use experienced log returns(in the paper)
tmur_2_full = expweightedavg(exret_full , I.nu); %2 to T_total, use experienced log excess returns(in case)

% Trailing 4-quarter sum of dividends
D_ttm_full = movsum(exp(d_full) , 4 , 'Endpoints' , 'discard');

PD_ttm_full = P_full(4 : end , :) ./ D_ttm_full; %For the last t periods
pd_ttm_full = log(PD_ttm_full);

tmud_full = expweightedavg(dg_full , I.nu);

%% 3. SIMULATED DATA

%% 3.1 Drop burn-in periods
% All data are re-labled from 1 to T as realizations. 
% Be careful with risk-free rate

%-------------------------------------
% Endowment Process
%-------------------------------------
cg = cg_full(end - T + 1 : end , :);
dg = dg_full(end - T + 1 : end , :);
c  =  c_full(end - T + 1 : end , :);
d  =  d_full(end - T + 1 : end , :);
cg_yoy = cg_yoy_full(end - T + 1 : end , :);
dg_yoy = dg_yoy_full(end - T + 1 : end , :);

%-------------------------------------
% State variable proxies
%-------------------------------------
tmu  =  tmu_full(end - T + 1 : end , :);
tmud = tmud_full(end - T + 1 : end , :);
tmur_1 = tmur_1_full(end - T + 1 : end , :);
tmur_2 = tmur_2_full(end - T + 1 : end , :);
tmur = tmur_1; % For now we use the experienced log returns specification

%-------------------------------------
% Price Quantity
%-------------------------------------
rf = rf_full(end - T + 1 : end , :);

P = P_full(end - T + 1 : end , :);
D_ttm = D_ttm_full(end - T + 1 : end , :);
pd_ttm = pd_ttm_full(end - T + 1 : end , :);
pd = pd_full(end - T + 1 : end , :);

ret = ret_full(end - T + 1 : end , :);
ret_ll = ret_ll_full(end - T + 1 : end , :);

ret_MS = ret_MS_full(end - T + 1 : end , :);

EXRET = EXRET_full(end - T + 1 : end , :);
EXRET_ll = EXRET_ll_full(end - T + 1 : end , :);

exret = exret_full(end - T + 1 : end , :);
exret_ll = exret_ll_full(end - T + 1 : end , :);

sub_ret = sub_ret_full(end - T + 1 : end , :); 

Sharpe = I.sharpeconst + I.sharpecoeff * tmu;

%% 3.2 Subjective growth expectation
t_fc = 20; % Number of periods to forecast
% Calculate coefficients
A_K = I.phi * (1 - I.nu + I.phi * I.nu).^(0 : t_fc - 1)';
A_K_s = I.lambda * A_K(2 : end) + (I.alpha - I.lambda) * A_K(1 : end - 1) - I.alpha;
A_K_sum = expweightedavg(A_K_s , I.alpha)/I.alpha;

% Calculate subjective expectation 
LTG = nan(T , N_sim); % Average subjective dividend growth expectations for the next 5 years
Y1G = nan(T , N_sim); % Average subjective dividend growth expectations for the next 1 year
sub_dg_temp = nan(t_fc , N_sim); % Store the next 20-quarter subjective dividend growth expectations at each point of time
for k = 1 : T
    % Separately calculate \tE_t\Delta d_{t+1}
    sub_dg_temp(1 , :) = I.lambda * (I.phi * tmu(k , :) + (1 - I.phi) * I.mu) - I.alpha * (d(k , :) - c(k , :) - I.mu_dc);
    
    % Calculate all the following periods
    sub_dg_temp(2 : end , :) = ...
    repmat(sub_dg_temp(1 , :) , t_fc - 1 , 1) .* repmat((1 - I.alpha).^(1 : t_fc - 1)' , 1 , N_sim) ...
  + repmat(tmu(k , :) , t_fc - 1 , 1) .* repmat(1 - (1 - I.alpha).^(1 : t_fc - 1)' , 1 , N_sim) ...
  + repmat(tmu(k , :) - I.mu , t_fc - 1 , 1) .* repmat(A_K_sum , 1 , N_sim);
            
    LTG(k , :) = mean(sub_dg_temp); 
    Y1G(k , :) = mean(sub_dg_temp(1 : 4 , :));
end
 
%% 4. STORE SIMULATIONS
SIMDATA = cell(length(Var_Names) , 1);
for i = 1 : length(Var_Names)
    SIMDATA{i} = eval(Var_Names{i});
end

end